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In soft matter systems the local displacement field can be accessed directly by video microscopy 
enabling one to compute local strain fields and hence the elastic moduli in these systems using a 
coarse-graining procedure. Here, we study this process in detail for a simple triangular lattice of 
particles connected by harmonic springs in two-dimensions. Coarse-graining local strains obtained 
from particle configurations in a Monte Carlo simulation generates non-trivial, non-local strain 
correlations (susceptibilities), which may be understood within a generalized, Landau type elastic 
Hamiltonian containing up to quartic terms in strain gradients (K. Franzrahe et. ai, Phys. Rev. E 
78, 026106 (2008)). In order to demonstrate the versatility of the analysis of these correlations and 
to make our calculations directly relevant for experiments on colloidal solids, we systematically study 
in detail various parameters such as the choice of statistical ensemble, presence of external pressure 
and boundary conditions. Crucially, we show that special care needs to be taken for an accurate 
application of our results to actual experiments, where the analyzed area is embedded within a 
larger system, to which it is mechanically coupled. Apart from the smooth, affine strain fields, the 
coarse-graining procedure also gives rise to a noise field (x) made up of non-affine displacements. 
Several properties of x ma y be rationalized for the harmonic solid using a simple "cell model" 
calculation. Furthermore the scaling behavior of the probability distribution of the noise field (x) 
is studied. We find that for any inverse temperature /3, spring constant /, density p and coarse- 
graining length A the probability distribution can be obtained from a master curve of the scaling 
variable X = x/3//pA 2 . 



PACS numbers: 62.20.D-, 82.70.Dd, 05.10.Ln 
I. INTRODUCTION 

Soft matter with its structural and elastic properties 
offers an attractive route to the design of new materials. 
In particular colloidal dispersions attract a lot of inter- 
est in this context. Surface chemistry or alterations in 
the composition of the solvent give an excellent control 
over the effective interactions in colloidal dispersions 
By definition colloids lie in the range of the visible spec- 
trum. Video microscopy |2| is therefore a straightforward 
means to gain information of the microscopic trajecto- 
ries of the components of the system under study. Thus 
microscopic, thermal (or Brownian) fluctuations can be 
resolved directly in real space, making colloidal disper- 
sions excellent model systems for the study of fundamen- 
tal questions of the statistical physics of soft condensed 
matter. Two dimensional colloidal dispersions, for ex- 
ample, have been used successfully in studies on melting 
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in two dimensions during the last decades [1, 0, H, @]- 
In this paper we focus on the mechanical properties of 
such systems and consider, in detail, how they may be 
obtained from the microscopic particle trajectories. 

Within linear elasticity, a solid in two-dimensions 
is described by eight unknown variables: the three 
stresses <7y, three strains and two components of 
the displacement field it,. Appropriately, there are also 
eight equations, namely the two equilibrium conditions 
dtjji/dxj + fi = (with fi, the forces per unit vol- 
ume within the body), three geometrical equations ey = 
{dui/dxj + duj/dxi) /2, and three constitutive equations 
Uij = Cijkitki- This set of equations may be solved for a 
given boundary condition in order to extract either the 
stresses or strains, given the elastic moduli, or the elastic 
moduli themselves, if the strains are known for a given 
stress configuration or vice versa. This manner of obtain- 
ing elastic moduli requires us to perturb the system using 
some external means e.g. laser tweezers Q- In contrast 
to this approach, one may calculate the tensor of elas- 
tic constants Cijki of a system from fluctuations of the 
microscopic strains obtained by a coarse-graining proce- 
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dure. Computing Cyjy in this way requires no external 
forces to be applied which may tend to change the very 
properties that are being measured @, ©, HI- Recently, 
this procedure was further extended in Ref. to obtain 
even the non-local elastic susceptibilities. 

The purpose of the present paper is twofold. Firstly, 
we present in great detail the analytic background of 
the coarse-graining procedure for obtaining elastic mod- 
uli and non-local susceptibilities (or strain-strain corre- 
lation functions) used in our earlier work [llj. In order 
to demonstrate the versatility inherent to the analysis 
of these non-local elastic susceptibilities, we present sys- 
tematic studies of a simple two dimensional lattice of 
particles connected by harmonic springs in the current 
paper. A comparison of the non-local susceptibilities in 
different statistical ensembles for various boundary con- 
ditions, system sizes and under different external condi- 
tions is given. Furthermore relations between the non- 
local susceptibilities and the elastic constants in systems 
surrounded by an embedding medium are derived. The 
proper interpretation of the correlation functions in such 
settings is discussed and visualized by use of the static 
susceptibility sum rule. Thus our aim in the present pa- 
per is to demonstrate various approaches in the analysis 
of the non-local strain correlation functions and to show 
how the analysis has to be adapted to the actual exper- 
imental situation. This study will thus greatly facilitate 
adoption of such techniques for routine analysis of exper- 
imental data at least for soft systems, which are close to 
being harmonic. 

Secondly, apart from the above stated intent to estab- 
lish a precise procedure for obtaining mechanical proper- 
ties from microscopic configuration data, we also aim, to 
study in some detail fundamental aspects of the coarse- 
graining procedure itself. For example, an immediate 
problem is the presence of particle configurations within 
the coarse-graining volume, which are not describable in 
terms of affine deformations of any reference lattice, e.g. 
incipient vacancy-interstitial pairs. This is true for all 
coarse-graining volumes larger than an unit cell. What is 
the effect of these configurations on elasticity and how do 
they influence mechanical behavior? Recently, there has 
been significant progress in the study of non-affineness in 
solid plasticity - especially in the context of rheological 
properties of amorphous materials and granular solids 
which show jamming behavior [13 . [l3l Il4j |. Localized 
non-affine regions consisting of particles capable of large 
reorientations have been shown to be involved in relax- 
ation processes in these systems. Is there an analog of 
such regions in an ideal, crystalline solid? A study of 
these fluctuations in ideal solids, as presented in section 
IIV1 may help us understand complex dynamics in solids 
better. 

The organization of this paper, together with a short 
summary of our main results is as follows. In section 
HIl we derive an analytic form of the non-local elastic re- 
sponse function, or compliance Xiji^i ^) (* = x iU)j which 
is defined as the strain £ij{f) produced at position r* due 



to a stress (Jij(r) at r. In order to do this, we consider 
a Landau expansion [TEj of the free energy in terms of 
the strains, keeping up to quartic terms in the gradients. 
Next in section IIIII we present Monte Carlo computer 
simulations of a harmonic crystal. The calculation of the 
local strain field corresponds to a coarse-graining proce- 
dure and allows us to construct the strain-strain correla- 
tion function, Gij which is related to the response func- 
tion via Xij = (ksT)~ 1 Gij. For a homogeneous solid 
without external load (i.e. (si(r)) = 0) the correlation 
functions are given by Gij(f*) — V < £i(0) £j{f) >■ The 
< ... > denote a thermal average (and in addition one 
over the choice of origin) and ksT is the Boltzmann con- 
stant times the temperature. We compare our results, 
obtained for a variety of ensembles and boundary condi- 
tions to that of the Landau theory. A common feature 
in experimental systems is the presence of an embedding 
medium, surrounding the analyzed region of the sample. 
The effects of such an embedding medium on the strain 
correlations are discussed and visualized by use of a sta- 
tistical sum rule. One of our significant results is that 
though the forms of the correlation functions and their 
limiting values as predicted by the Landau theory are re- 
produced, the Gij obtained from simulations through our 
coarse-graining procedure differ by an additional back- 
ground contribution which is not negligible. In section 
IIVI we argue that this is a consequence of non-affine dis- 
placements, which are not considered in the ansatz for 
the Landau theory. The amount of non-affinity in a 
given configuration can be quantified by calculating x, 
the deviation of the actual configuration from one ob- 
tained from an affine transformation of the reference lat- 
tice. We analyze this non-affineness in detail and show 
that the probability distribution of the non-affine field 
P(x) can be computed within a simple "cell model" ap- 
proximation. P(x) shows well defined scaling properties 
with the spring stiffness / and the coarse-graining length 
A. The auto-correlation function for x is shown to be 
short ranged decaying rapidly for distances much larger 
than the coarse-graining length. Finally, we conclude our 
paper indicating future directions for research. 

II. LANDAU THEORY FOR THE STRAINS 

The two dimensional elastic continuum described by 
the linear elastic Hamiltonian, 

PH = 1 -JdfC ijkl e ij e kl (1) 

is perpetually in a critical state [l5l |. The displacement 
correlations (u(r)-u(f')) decay algebraically and the solid 
shows quasi-long ranged order with the elastic suscepti- 
bilities diverging logarithmically with system size L. For 
all practical purposes, however, this weak divergence may 
be ignored and non-zero elastic moduli may be defined 
and computed. In real solids, an upper length scale cutoff 
is set by the typical distance between defect pairs. 
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The fact that the solid state is critical also implies that 
the Hamiltonian Hq in Eq.|T|) is a fixed point Hamilto- 
nian which should be invariant under a coarse- graining 
procedure, unless topological defects such as dislocations 
are present. This, again, is not strictly true, as we shall 
demonstrate in section iHfl This is because in a molecu- 
lar system, any Hamiltonian such as TIq is realized only 
in an approximate, discretized sense, the displacements 
u being carried by the individual particles. Anticipating 
some of our results in section IIII] we use the following 
dimensionless Landau functional Il5l. fl6j : 



fiT = 



1 



d 2 x 



3 

E 

i=l 



{a l e l 2 + c l (Ve i ) 2 + ^(V 2 e l ) 2 }(2) 



Here with (i = 1 — 3) the dimensionless constants a% 
are the elastic moduli of the system, while a and c\ are 
phenomenological coefficients and the three strains are 
given by: ei = e xx + e yy , describing pure volume changes; 
62 = ^xx — tyy, describing deviatoric shear strains and 
e 3 = \{t X y + describing pure shear strains. The 
phenomenological coefficients ci have the dimension of a 
length 2 . Thus we interpret £ e j j ~ Jcl as a correlation 
length. Note that unless noted otherwise throughout the 
paper all lengths are given in units of a, the lattice param- 
eter of the underlying triangular lattice in the simulated 
systems. 

The terms quadratic in the strains e\ represent the 
local part in this ansatz. Non-local contributions are 
included via the gradient terms. Note that the Lan- 
dau functional Eq. (J2j) should be strictly valid for exci- 
tations of wavelengths longer than a short-wavelength 
cutoff. Short wavelen gth excitations are suppressed by 
the gradient terms [IRTl6j and we have ignored the pos- 
sibility of defects. 

For non- uniform strains, spatial fluctuations of the 
strains couple and only one of the strain variables ej in 
this ansatz is an independent variable as we show below. 
Firstly, all forces within and on the system must can- 
cel for the system to be in thermodynamic equilibrium, 
for a solid under zero external stress; the 
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dxj 
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case of external hydrostatic stress is presented in section 
IIII CI The stress tensor can be obtained directly from 
Eq.®, as CTjj = S g ^' F . For a two-dimensional crystal this 
condition reads in Fourier space: 

k x a x e x + k x a 2 e 2 + k y a 3 e 3 = (3) 
k y a x e x - k y a 2 e 2 + k x a 3 e 3 = (4) 

In addition St Venant's compatibility condition must 
be considered in the calculations, which ensures an 
unique relation between the displacement field u and the 
strain fields e,j. For a two-dimensional crystal this sim- 
plifies to 



fc 2 ei + (fc 2 - kl)e 2 + ik x k y e 3 = 



(5) 



in Fourier space. With the help of these three conditions 
the kernels Qij relating the strain variables = Qijfij 



can be derived. The resulting relations are given in detail 
below: 



e2 = 



ei 



4ai + 2Q3 ^ ( k * k v ) g 3 = Q 23 ~e 3 
ai+a 2 ) \k 2 x -k 2 y ) 

2a 3 - 4a 2 \ / k x k y \ _ ~ _ 
e 3 = <yi 3 e 3 



ei 



ai + a 2 J \ k 2 
as - 2a 2 \ / kl - k 2 y 
lax +03/1 k 2 



62 = Qi2e 2 



Note that the properties of the correlation functions are 
set by the wave vector dependence of the kernels Qij . 

We shall also need kernels relating the strains to the 
local, microscopic rotations given by the anti-symmetric 
part of the strain tensor 9 = (du y /dx — du x /dy) /2. Lo- 
cal rotations 9 are related to the deviatoric strains e 2 
and the pure shear strains e 3 . These strains are cou- 
pled in case of the presence of a surrounding, embedding 
medium. Starting from the definitions of the strain vari- 
ables the partial derivative d 3 9/dx 2 dy can be expressed 
solely in terms of e 2 , e 3 and 9. Thus one obtains the fol- 
lowing relation between the three shear strain variables 
in Fourier space: 



9{k 2 x + fc 2 ) = (fc 2 



fc 2 )e 3 - k x k y e 2 



We may now derive a new kernel relating 9 to e.g. the 
pure shear strain e 3 using the previously derived kernel 
Q 23 . This results in: 



k, r 
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fU2 _ l,2\2 i £.21,2 ( 4ai+2a 3 A 
[K x K y ) -t- K x K y y Qi+02 j 



= Q 



As the considerations in section ITlID I will show, it will be 
helpful to define the strain variable e 2 g = 29 and analyze 
its two-point correlation function in embedded systems, 
such as the colloidal crystal discussed in 

Switching to a discretized notation, appropriate for 
comparison with our simulations, so that e\(x) — > e\ ^ 
and x — ► i?^, where m (a tuple of integer lattice indices) 
specifies the position on a coarse-graining, square mesh 
the free energy can be rewritten as, 



Here v z = V z /a 2 is the dimensionless volume of the 
coarse-graining cell and the summation runs over all N 
cells. One may use the above relations to finally express 
(3F as a harmonic functional of only one of the strain 
components. This allows the direct calculation of the an- 
alytic form of the two-point correlation functions, as the 
partition function factorizes (e.g. [I3|)- Here we choose 
as an example the strains e\ and obtain the following 
expression for the functional, written as a sum over the 
wave vectors: 
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fiT = 



5?<« 



k 2 Cl + fc 4 ci + (a 2 + c 2 k 2 + c' 2 k 4 )(Q 21 (k)) 2 + (a 3 + c 3 k 2 + c 3 fc 4 )(Q 31 (fc)) 2 }^ - 



r 



with v = V/a 2 , the dimensionless volume of the system. 
It is now straight forward to calculate the two-point cor- 
relation function starting from its general definition and 
taking into account the fact, that the average strains 
(ej(f)) in a crystal, which is under no load, is zero. 



teWei(0>-<ei(f)>te(0) 



E 



A-A-' 



A(k-f)(k'-r))l~~^ 



The calculation yields the following relation for the aver- 
age of the Fourier coefficients: 

Where, A{k) = a a + k 2 c x + k 4 c[ + (a 2 + c 2 k 2 + 

4fc 4 )(Q 21 (fc)) 2 + (a 3 +c 3 fc 2 + C ^ 4 )(Q 3 i(fc)) 2 . AsG u (k) = 
(|eg| 2 ) /v one can now write down the analytic form of the 
strain-strain correlation functions or rather their inverse 
in detail. For i — 1,2,3 the structure of the strain-strain 
correlation functions is the same, while for i = 29 it dif- 
fers slightly: 

G^fe^O)" 1 = a i + fc 2 c J + fc 4 ^ (6a) 

3 

+ Y, (a ] +c 3 k 2 +c' ] k i )(Q ]1 (k)) 2 
G^Or 1 = a z 



G 
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(k^O)- 1 = (aa + c 3 k 2 + c' 3 k 4 



(6b) 



G 2 g 2 g(0) 1 



+ X>j + cjk 2 + c> 4 )(Q i3 ) a ) (Q 3 20)' 

3=1 
«3 
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A. Properties of the strain-strain correlation 
functions 

The analytic strain-strain correlation functions Gu(k) 
(the inverse of which were given in Eq.©) are plotted 
in figure [TJ The set of parameters Of, Cj and c£ used in 
figure [T] were obtained from Monte-Carlo simulations in 
the NVT ensemble with periodic boundary conditions of 
a harmonic triangular lattice, to be discussed in section 



1111 Bl While the deviatoric and shear strain correlation 
functions G 22 and G 33 have four-fold symmetries, the 
correlation function of the dilatation Gn has an eight- 
fold symmetry. The correlation functions may be inter- 
preted as the response of the system to a localized per- 
turbation at the origin. This perturbation is either a 
dilatation, a deviatoric shear or a pure shear. The de- 
formation of the solid may be decomposed as a super- 
position of the eigenmodes of the system. The eigen- 
modes for a square box, are plane waves with polariza- 
tions either longitudinal or transverse to the coordinate 
axes with the eigenfrequencies forming a discrete spec- 
trum: u> nm = ^jjCa^n 2 + m 2 with n, m G No- Thus the 
wave vector of the eigenmodes along the diagonal, i.e. 
(n = to), exhibits a four- fold degeneracy, while those par- 
allel to the coordinate axis, i.e. n ^ m and n/m^O, 
have an eight-fold degeneracy. A local dilatation as per- 
turbation results in a superposition of eigenmodes with 
four- as well as eight-fold degeneracy. This leads to the 
eight-fold rotational symmetry visible for the strain cor- 
relation function Gn. In contrast to this the two possi- 
ble shear perturbations will excite elastic waves that are 
superpositions of exclusively eigenmodes with four-fold 
degeneracy. For this reason the corresponding correla- 
tion functions G 22 and G 33 exhibit only a four-fold rota- 
tional symmetry. As was discussed in 18j, the presence 
of defects, breaks the rotational symmetries of the strain 
correlation functions. 

The details of the structure of the correlation functions 
arc dominated by the dependence of the kernels Qij on 
the wave vector k, especially the cases k x = k y , k x — ► 
while k y ^ and k y — > while k x ^= 0. In particular, we 
obtain the following relations: 
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k x ^0 
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k x ^0 



For the kernel Q 3 2 g relating e 3 to e 2 g the behavior along 
the specific directions in Fourier space for k — > can be 
extracted from the behavior of the correlation function 
G 33 (k) and G 22 (k). Along the coordinate axis the kernel 
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FIG. 2: (Color online) Cuts of the analytic strain-strain cor- 
relation functions along specific directions in Fourier space. 
The set of parameters used for plotting is: a\ = 100, a 2 = 50, 
a 3 = 200 and ci = 54.3, c 2 = 37.4, c 3 = 118.2, c[ = -86.1, 
c 2 = —1.3, c 3 = —18.8. The correlations lengths were ob- 
tained from a Monte Carlo simulation of a harmonic trian- 
gular lattice in the NVT ensemble with periodic boundary 
conditions. Gn(fc) is shown along the direction k y = 2k x , 
C?22(fc) along the diagonal k x — k y , G 33 (fc) along the fc H -axis 
and G202e(k) along the diagonal k x = k y and along the k y - 
axis. 



shown to equal 1. Upon insertion of these relations into 
the equations for the inverse of the correlation functions 
their behavior for these limiting cases can be extracted: 
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FIG. 1: (Color online) The analytic form of the strain-strain 
correlation functions. For each function a surface plot and 
next to it a density plot are shown. In the density plot the 
maxima are white, minima are black. The set of parameters 
used for plotting was obtained from a Monte Carlo simulation 
of a harmonic triangular lattice in the NVT ensemble with 
periodic boundary conditions: a\ — 97.4, at — 48.1, a 3 = 
198.2, ci = 54.3, c 2 = 37.4, c 3 = 118.2, ci = -86.1, c' 2 = -1.3 
and c 3 = -18.8. a) Gn(fc), b) G 22 (fc) , c) G 33 (fc) and d) 

G2S2fl(fc). 
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relating e 3 to §20 becomes a constant 
Q3 29 = 



-1/2 
1/2 



for k 
for k 



u 



ky^O 

k x ^0 



Thus the continuous behavior of the correlation function 
G 33 (k) for k — ► along the coordinate axis carries over 
to the correlation function G282s(k). The behavior along 
the diagonals k x = k y can be extracted from the behavior 
of the product of the kernels Q\ 3 Q\ 2 e ■> which can be 
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These considerations show, that in certain directions in 
Fourier space the shear strain variables become indepen- 
dent from each other and are continuous for k — > . 
These are the directions, along which a fit will give direct 
access to the elastic constants and correlation lengths of 
the system. Figure [5] shows cuts in Fourier space which 
correspond to these specific directions for the strain cor- 
relation functions Gn{k). For the correlation function of 
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the pure shear strain G 22 (k) a cut along the diagonals 
is shown, while for the deviatoric strain correlation func- 
tion G3s(k) a cut along the coordinate axis is shown. The 
strain correlation function for the dilatation Gn(fc) in 
contrast is not continuous for k — > 0. If one considers for 
example the direction k y — 2k x the kernels Qij relating 
the strain variables turn into constant weighting factors: 

= ~ (ct) (f) and <33i = (iSl) (I)- Thus 
along the direction k y = 2k x the inverse of the strain 
correlation function for the dilatation G^j (fc) for k — > 
is given by 



G-^fc^O) - [ ai +a 2 



( 2oi + 03 \ / 5 



«3 



\a 3 — 2a 2 J \3 y 
2 s 



ai + a 2 
4a 2 ~2a 3 J \2 



G^(k = 0) = ai 

So the strain correlation function for the dilatation 
Gn(fc) exhibits a pronounced discontinuity for A; — > 0. 
In order to illustrate this fact, consider the set of param- 
eters used in the simulations of a harmonic triangular 
lattice (to be discussed in section llllj) . The choice of the 
spring constant sets the elastic constants of the system 
under consideration to a\ = 100, a 2 — 50 and 03 = 200. 
For k — > one has G^(k — ► Q) « 3025, which is ap- 
proximately 30 times as much as the value for k = 0, i.e. 
Gf x (fc = 0) = ai = 100, set by the bulk modulus - 
showing that non-uniform dilations tend to be severely 
penalized in this solid. 

Nevertheless, provided that the value of the bulk mod- 
ulus is determined e.g. from Gu(fc = 0) the coefficients 
ci and c[ can also be obtained by fitting one of the cor- 
relation functions along a cut in Fourier space. So in 
principle all 9 parameters of the free energy functional 
can be determined from an analysis of the strain-strain 
correlation functions. Like the correlation function of the 
dilatation Gn(fc) the correlation function of the micro- 
scopic rotations shows an eight- fold rotational symmetry. 
Unlike Gn(k), however, G 2 e 2 e(k) is continuous for k — > 
along the coordinate axis and the diagonal (compare fig- 
ure H]) . Therefore fits along these directions can be used 
to determine the elastic constants a 2 and 0,3 as well as 
the coefficients c 2 , c' 2 and C3, c 3 . 

We shall next discuss the result of a coarse-graining 
procedure, which attempts to obtain these correlation 
functions and therefore the parameters of the Landau 
free energy functional from Monte Carlo simulations of 
the harmonic lattice. The coefficients of all the second 
and fourth order terms involving gradients of strain are 
found to be non-vanishing showing that coarse-graining 
generates these higher order non-local terms in the free 
energy. 



III. MONTE CARLO SIMULATIONS OF A 
HARMONIC CRYSTAL 

The analysis of a harmonic crystal is convenient for a 
comparison with the results of the Landau theory pre- 
sented in the last section, since the elastic moduli can be 
directly calculated from the spring constants. We con- 
sider a harmonic triangular lattice with a Hamiltonian 

H = k B T{f/2)YZ,n=iiVm ~ r n \ - a) 2 where / is the 
spring constant and a the lattice parameter of the trian- 
gular lattice. The elastic moduli arc related to the spring 
constant / via: K = a\ = (V3/2)/, /j, = a 2 = (V3/4)/ 
and 4/x = 03 = >/3f. Furthermore the harmonic trian- 
gular lattice has been shown to be a successful model 
for the interpretation of experiments on colloidal crys- 
tals [ill] . It is modeled by N point-particles each of 
them hard-wired by spring constants / to the six nearest 
neighbors. We have carried out Monte Carlo simulations 
in the constant NpT and NVT ensembles with periodic 
boundary conditions. We also mention briefly results for 
a system with open boundary conditions which were pre- 
sented elsewhere [Tl|] . Next the influence of hydrostatic 
pressure is analyzed by Monte Carlo simulations in the 
constant NpT ensembles with periodic boundary condi- 
tions. Finally we consider the effect of a surrounding 
elastic medium and finite size effects in order to make 
contact with experiments on colloids. 

The knowledge of the configurations and the reference 
lattice allows for a direct calculation of the displacement 
field u(r) . In order to calculate the corresponding strain 
field partial differentials of the displacement field have 
to be calculated. We follow the procedure by Falk and 
Langer [l2j and calculate the strain field by minimizing 
the error in the affine transformation that relates the 
actual configuration {r} to the reference lattice {R}. 

r = R + u(R) = (l + e)R 

The mean-squared error in this mapping x is thus a mea- 
sure of how well the given situation can be described 
within the framework of linear elasticity theory and quan- 
tifies the non-affinity of the given displacement field. Falk 
and Langer [12| analyzed the temporal development of 
strains. Here we use an analogous definition for \ in 
thermodynamics equilibrium, evaluating the strains and 
non-affineness with respect to the reference lattice: 



N B 



Xiro) = E E 



E(% + e«)[<-4] 



Here rb is the position, at which the strains are to be cal- 
culated, and Nb is the number of neighboring particles 
considered. This corresponds to a coarse-graining proce- 
dure, in which Nb is set by the choice of coarse-graining 
length A, i.e. cutoff radius within which particles are 
considered in the calculation. For the results presented 
in this section, we have used a cutoff radius of A = 1.3 
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ai 


a 2 




calculated from / 


100 1 pa 


50 /(3a 2 


200 1 pa 


from fluctuations 

Ul 11 


QQ Q 1 ftn^ 
ao.ij I fja 


AQ A 1 a n 2 
'ry.'r 1 fja 


1 Qfi 7 /fln 2 
iyu. / / pa 


from G(k = 0) 


96.8 /Pa 2 


48.6 /Pa 2 


190.8 //3a 2 


from fits of G{k 0) 




49.1 /Pa 2 


195.7 /Pa 2 



TABLE I: A comparison of the elastic constants as calculated 
from the spring constant / and as obtained by use of various 
methods from data of a Monte Carlo simulation of a har- 
monic triangular lattice with N = 3120 particles in the NpT 
ensemble at zero hydrostatic pressure Pa 2 p = 0.0 and spring 
constant Pa 2 f = 200/^3- 



resulting in Nb = 6. In section HVl we present some sys- 
tematics showing how some of our results depend on the 
coarse-graining length A. 

In the calculation of the strain-strain correlation func- 
tions a second coarse-graining step is employed, when 
mapping the triangular lattice to a square mesh. This 
facilitates the numerical Fourier transformation of the 
calculated real space correlation functions. The wave 
vectors are limited to the first Brillouin zone, i.e. kj € 
[— f~if~\t with j = x,y. Here l m represents the lat- 
tice parameter of the coarse-grained, square mesh. Care 
must be taken in the choice of l m to keep the coarse- 
graining volume large enough so that artifacts due to the 
discreteness of the triangular lattice (and insufficient av- 
eraging) are avoided. In most of the results presented 
here we used Z m = 2.25a, where a is the lattice param- 
eter of the original, triangular lattice. Lastly, one also 
needs to be careful about correcting for global rotations 
and translations of the lattice especially for the case of 
open boundary conditions so as not to introduce artificial 
sources of error. 

Simulations of the harmonic triangular lattice were 
done for three system sizes N = 3120, 4736 and 5822. We 
first discuss the results for simulations with spring con- 
stant f3a 2 f — 200/\/3 in the NpT ensemble with periodic 
boundary conditions (and external pressure p = 0) and 
the NVT ensemble with periodic boundary conditions. A 
crystal with open boundary conditions was discussed in 
detail in [TTI ] and will only be mentioned briefly. In what 
follows we analyze the influence of a hydrostatic, external 
pressure and of a surrounding, embedding medium. 



NpT ensemble with periodic boundary 
conditions 





FIG. 3: (Color online) Strain-strain correlation functions of 
a harmonic triangular lattice at zero hydrostatic pressure as 
obtained from Monte Carlo simulations in the NpT ensemble 
with periodic boundary conditions. Results for a system with 
N — 3120 particles and a spring constant Pa 2 f = 200/\/3 
are shown. For each function a surface plot and next to it a 
density plot are displayed. In the density plot the maxima 
are white, minima are black. 



semble, is saved in the transformation matrix h. One of 
the advantages of this implementation is, that from the 
fluctuations of the transformation matrix h the fluctua- 
tions of the strain tensor can be calculated directly. The 
strain tensor is related to the transformation matrix via 



For the simulations of the harmonic crystal in the NpT 
ensemble we use the algorithm of Parinello and Rahman 
[20| . Here the information on the actual shape of the 
simulation volume, which is free to fluctuate in this en- 



where h is the transformation matrix of the reference 
lattice and G = h T h contains the information of the ac- 
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FIG. 4: (Color online) Cuts of the strain-strain correlation 
functions in the NpT ensemble with periodic boundary con- 
ditions at zero hydrostatic pressure along specific directions 
in Fourier space, a) Gu(k) along k y = 2k x , b) G?22(fc) along 
k x = k y , c) G33(fe) and d) G262e{k) along the fc^-axis. For 
a comparison the results for different system sizes are dis- 
played. The horizontal broken lines mark the value expected 
from theory for the correlation functions at k = 0. 





ai 


ai 


a?, 


calculated from / 


i nn / ft ~ 2 

1UU 1 pa 


rn / ft n 2 

ou / pa 


onn / ft ~ 2 

2UU 1 pa 


using Squire et al. [211 


97.4 /Pa 2 


48.3 /Pa 2 


194.8 /Pa 2 


from fits of G(k / 0) 




48.1 /Pa 2 


198.2 /Pa 2 



TABLE II: A comparison of the elastic constants as calcu- 
lated from the spring constant / and as obtained by use of 
various methods from data of a Monte Carlo simulation of a 
harmonic triangular lattice in the NVT ensemble with peri- 
odic boundary conditions with N = 3120 particles and spring 
constant Pa 2 ) = 200/^3. 



shear modulus as it is given in table U and the coefficients 
c 2 = 34.3, c 2 = -0.6, c 3 = 114.1 and c' 3 = -17.6. So the 



elastic correlation lengths £ e ;^ ~ Jci are approximately 
6 and 11 lattice parameters respectively. 

Figure |4] d) shows cuts along the coordinate axis of 
Gie 2${k). In these simulations the system as a whole is 
not an embedded system and is not free to rotate. Thus 
we cannot obtain the elastic modulus directly from the 
value of G 2 g 28 at the origin. This situation is different 
in a solid, which is embedded in a larger volume, as will 
be discussed in section ITlIDl 



tual shape of the simulation volume. Table U shows a 
comparison of the elastic moduli for the harmonic, trian- 
gular lattice as they are expected for the chosen spring 
constant / and the values as they are obtained from the 
simulations in the NpT ensemble by analyzing the fluc- 
tuations of the simulation volume. 

Figure [3] shows the strain-strain correlation functions 
in Fourier space as they are obtained from the simula- 
tions in the NpT ensemble with periodic boundary con- 
ditions. The eight-fold rotational symmetry in Gn(fc) 
is not resolved. The shear strain correlation functions 
show clearly a four-fold rotational symmetry, as was ex- 
pected from the analytic predictions. Cuts along various 
directions in Fourier space of these functions are plotted 
in figure Q] for the three system sizes. As these cuts in 
Fourier space show, there is no systematic dependence 
on the system size in the correlation functions. 

The discontinuities in the correlation functions are vis- 
ible in figure [3] and |U Nevertheless the extreme disconti- 
nuity one expects to observe in Gn(fc) from the analytic 
predictions is reduced to a factor of approximately 1.5 
instead of 30, as a comparison of the cuts in figure [2] and 
in figured] a) shows. This indicates that there might be 
excitations in the system that are not captured by the 
assumption of purely affine strains. Along the cuts, for 
which G 2 2(fc) and G33(k) are continuous for k — > 0, fit- 
ting with a generalize Lorentzian profile yields the elastic 
constants and via the coefficients c, the elastic correlation 
lengths. For the system with N — 3120 we obtain the 



B. NVT ensemble with periodic boundary 
conditions 

Below, we describe simulations in the NVT ensemble 
with periodic boundary conditions at a reduced density 
of g* — 1.0. Table HT1 lists the elastic constants calcu- 
lated with the fluctuation method given by Squire et. 
al. [2l[. These authors also give a formula for calculat- 
ing the stress tensor. An evaluation of the data yields 
Pxy = ®yx — 0.0 and from the trace of the stress ten- 
sor j3a 2 p — — \{p xx + <Tyy) w 0.1. So we verify that the 
simulations represent a solid at approximately zero hy- 
drostatic pressure. 

In this ensemble the k = values of the correlation 
functions cannot be used to calculate the elastic con- 
stants directly. We are simulating an undeformed state 
of the crystal, thus the integral over the fluctuations of 
the strains over the complete simulation volume tends 
to zero in this ensemble. Therefore only fits along the 
directions, for which the correlation functions are contin- 
uous for k — * 0, give access to the elastic constants in 
this ensemble. As Gn(fc) has no such direction, the bulk 
modulus cannot be obtained in this way. From these 
considerations one expects the strain-strain correlation 
functions in the NVT ensemble to differ from those in 
the NpT ensemble for small absolute values of k. Figure[5] 
shows surface plots and density plots of the strain-strain 
correlation functions in Fourier space. As in section lHI Al 
the anisotropies are recovered well, except that the eight- 
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FIG. 5: (Color online) Strain-strain correlation functions of 
a harmonic triangular lattice as obtained from Monte Carlo 
simulations in the NVT ensemble with periodic boundary 
conditions. Results for s system with N = 3120 particles and 
a /3a 2 f = 200/^3 are shown. For each function a surface plot 
and next to it a density plot are displayed. In the density plot 
the maxima are white, minima are black. 
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FIG. 6: (Color online) Cuts of the strain-strain correlation 
functions in the NVT ensemble with periodic boundary con- 
ditions along specific directions in Fourier space, a) Gn(fc) 
along k y = 2k x , b) G22(k) along k x — k v , c) G32,(k) and d) 
G292e{k) along the fc^-axis. For a comparison the results for 
different system sizes are displayed. The horizontal broken 
lines mark the value expected from theory for the correlation 
functions at k — 0. 



from the coefficients: c 2 — 37.4, c' 2 — —1.3, C3 = 118.2 
and c 3 = 18.8. So consistent to the results obtained from 
the simulations in the NpT ensemble £ e j^ ~ 6 and 11 
lattice parameters respectively. Fitting e.g. G 2 2(k 7^ 0) 
along the direction k y = 2k x allows the determination of 
the coefficients c\ — 54.3 and c[ — —86.1. Thus correla- 
tions of volume fluctuations decay over approximately 7 
lattice parameters. 

The harmonic, triangular lattice in the NVT ensemble 
was also analyzed with open boundary conditions. The 
results were discussed in detail in 111!. 



C. The influence of hydrostatic pressure 



fold rotational symmetry of Gu(k) is not resolved. The 
expected discontinuous jump to Ga(k = 0) = (i = 1, 
2, 3) is clearly visible. Besides this, the correlation func- 
tions coincide with those obtained in the NpT ensemble, 
as one can see by comparing figure [4] and [6l showing 
the same cuts in Fourier space for the various correlation 
functions. The elastic constants a 2 and 03, as they are 
obtained from fitting the strain-strain correlation func- 
tions, are listed in table [TTJ They fall within 3 — 4% of the 
theoretical values and have thus the same accuracy as the 
values obtained via Squire's fluctuation formulae [2lj . In 
addition the elastic correlation lengths could be obtained 



How does an external, hydrostatic pressure - i.e. a xy = 
o~yx = and a xx — <7 yy = —p - influence the strain-strain 
correlation functions? Simulations of a harmonic, trian- 
gular lattice with spring constant (3a 2 f = 200/v / 3 sub- 
jected to an external, hydrostatic pressure 13a 2 p — 20/ V3 
show, that the shape of the correlation functions is not 
affected. The NpT ensemble was chosen for this study. 
Strains were calculated with respect to the average lat- 
tice positions, that is the compressed lattice. The lat- 
tice parameter of this reference lattice a' — 1.010465 is 
smaller than the lattice parameter a — (2/v / 3) 1 ^ 2 in the 
zero-pressure simulations. Therefore for a comparison 
with the theoretical values, which were given in units of 
(3a 2 , the a* as they are obtained e.g. from Gu(k = 0) 
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FIG. 7: (Color online) Cuts of the strain-strain correlation 
functions in the NpT ensemble with periodic boundary condi- 
tions along specific directions in Fourier space. Shown are the 
correlation functions obtained in a system without (crosses) 
external pressure and a system with f3a 2 p = 20/\/3 (cir- 
cles). Both systems have N = 3120 and a spring constant 
Pa 2 f = 200/^. a) Gn(fc) along k y = 2k x , b) G 22 (fc) along 
k x = k y , c) G3s(k) and d) G2B20 (k) along the fc H -axis. 



must be rescaled to these units. Simulations were run 
for a system with N = 3120 particles. For a direct com- 
parison of the correlation functions in systems with and 
without a hydrostatic pressure cuts in Fourier space of 
the Gaik) (i — 1,2,3,20) are shown in figure [71 These 
show clearly the shift in the absolute values and the per- 



sistence of their shape. When relating the parameter eij 
to the elastic moduli of the system, one has to recall, 
that the strain-strain correlations are related to the stiff- 
ness tensor Bijki , which is defined via the stress-strain 
relations. These relate the variation of stress to the 
variation of strain, to first order in the strains, for the 
case, that an arbitrary initial configuration {R} is trans- 
formed to a final configuration {r} by an applied uni- 
form stress. Thus B ijM ({R}) = (dajj({r}) / 'de k i) { jh = 

\ {&il5jk + Vjl&ik + VikSjl + VjkSil — 2&ijSkl) + Cijkl [22|. 

The stiffness tensor explicitly depends on the applied 
stress. Only for the case that cry = is it equivalent 
to the tensor of the elastic constants Cyjy . For the cal- 
culation of the elastic moduli the following combinations 
are needed: 



2 (B X XXX ~\~ B XX yy^ 



(BxXXX B XX yy^ 



B, 



1 



1 

2 

K 

77 {Cxxxx C X xyy 

) - o (p + p) 

(j>-p 

Cxyxy V M P 



Thus the bulk modulus can be directly obtained from 
Oi, while from G*22(fc) one extracts 02 = /i — p and from 
G*33(fc) one obtains 03 = 4(fi — p). The values of the 
elastic moduli calculated according to this scheme are 
given in table Hill They lie within 4.4% of the theoretical 
values. 



D. Effects of an embedding medium and finite-size 



calculated 
from / 


100 /Pa 2 


H — £12 + P 

= 50 /Pa 2 


4/i = a 3 + 4p 
= 200 /Pa 2 


from 
G(k = 0) 


88.8 /Pa' 2 


a 2 = 36.0 /Pa' 2 


a 3 = 137.6 /Pa' 2 


rescaled 
values 


100.4 /Pa 2 


a 2 = 40.7 /Pa 2 
-> ii = 52.2 /Pa 2 


a s = 155.6 /Pa 2 
-> 4/i = 201.8 /Pa 2 



TABLE III: The elastic moduli of the harmonic system calcu- 
lated from the spring constant /3a 2 / = 200/\/3 in comparison 
to simulation results. Listed are the elastic constants as ob- 
tained from Ga(k = 0) for a simulation in the NpT ensemble 
with periodic boundary conditions of a harmonic, triangular 
lattice with Pa 2 p — 20/\/3- For comparison with the theoret- 
ical values the at obtained in units of the lattice parameter a' 
of the compressed reference lattice have to be rescaled to the 
lattice parameter a of the zero-pressure reference lattice and 
the relation between the stiffness tensor Bijki and the tensor 
of elastic constants djki in a system with Oij ^ need to be 
considered. 



Often the strain-strain correlations cannot be evalu- 
ated over the complete crystal. In experiments as for 
example a two dimensional colloidal crystal 0, [ll| con- 
figurational data is taken via video microscopy. Here the 
area accessible to the video camera is far smaller than the 
complete sample size. In these cases only a sub-system 
embedded in a larger continuum is analyzed. As Zahn 
et al. Q noted the presence of an infinite, embedding 
medium alters the relation of the strain fluctuations to 
the elastic moduli. 

The strain-strain correlation functions are the response 
functions to a strain perturbation at the origin. Figure [8] 
shows schematically the resulting displacement field and 
strain fields for the cases that this perturbation is a) a 
dilatation and b) a rotation. The connection between the 
strain correlations and the elastic moduli was derived un- 
der the assumption, that the considered functional of the 
free energy accounts for the free energy of the complete 
system (equipartition theorem). For the case that the 
volume Vb over which the strain-strain correlation func- 
tion are calculated is not equal to the complete system 
volume V this assumption is not fulfilled any more. As 
can be seen in the schematic plots of the strain fields in 
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FIG. 8: (Color online) A schematic drawing of the displace- 
ment field u(f) as it results from a given perturbation at the 
origin: a) a dilatation and b) a rotation a the disk at the 
origin. If the analyzed volume Vb (broken red line) does not 
coincide with the volume of the complete system V (black 
line), the energy needed for the displacements in the embed- 
ding medium (V — Vb) must be taken into account in the 
interpretation of the correlation functions. 



figure [5] the energy related to the resulting strain field 
outside Vb cannot be neglected for Vb ^ V. 

Following the argument by Zahn et al. j9j], but con- 
sidering a finite embedding continuum, we show that the 
influence of the surrounding medium on the strain fluc- 
tuations within the analyzed volume Vb depends on the 
relative size of Vb in comparison to the complete system 
volume V, i.e. the ratio of Vb /V. For the derivation of 
the formulae we consider first a homogeneous dilatation 
of a disk Vb = ^R\ m a surrounding medium of volume 
V = irR 2 and second a pure shear, which can be realized 



by a rotation by an angle of the disk with volume Vb ■ 
For these considerations we work in polar coordinates, 



where we have e± = (e r 



£ w ), e 2 = (e rr — e vv ) and 
£3 = trip- In both cases considered here it is assumed 
that the displacement field on the boundary of the com- 
plete system is given by u(f = Rboundary) = 0. 



1. Homogeneous dilatation 

An isotropic expansion of a disk embedded in a finite 
medium is given by: Rb — > Rb + Ar. The resulting 
displacement field in polar coordinates is given by: 



= 




for r < Rb 
for r > Rb 



From this it is straight forward to calculate the re- 
sulting strain field and consequently the Free Energy 
density / = \ [K(e rr + e w ) 2 + fi ({e rr - e w ) 2 + 4-e 2 .^)] 
of the system under load. Thus the total energy 
needed for such an expansion in a finite system of vol- 
ume V = ttR 2 is given by E = \^ ^ rdipdr f = 

{V B /2){AV B /V B f[K + ii{l-{V B /V))}. For such a 
system equipartition tells us thus, that the strain fluc- 
tuations are no longer set by the bulk modulus of the 
system, but aquire in the embedded system a term de- 
pendent on the shear modulus and on the ratio V B /V: 

knT , „, 1 



Vb 



Therefore the strain-strain correlation function Gn(r) 
no longer provides access to the bulk modulus, but to 
a V B /^-dependent combination of bulk and shear mod- 
ulus. 



2. Pure shear 

A rotation of the disk as a rigid body within the em- 
bedding medium by an infinitesimal angle Aip changes a 
given orientation 99 to ip + Aip. The resulting displace- 
ment field is given by: 









Ap- 



for r < Rb 
for r > Rb 



From the corresponding strain field the Free Energy den- 
sity can be determined and integration over the com- 
plete system yields the energy required for such a rota- 
tion: E = n(2Aip) 2 V B (1 - (V B /V)) /2. In case of in- 
finitesimal rotation angles Aip this angle can be identi- 
fied with the anti-symmetric part of the strain tensor 
6 = (duy/dx — du x /dy) /2. Equipartition relates the 
fluctuations in eie to the shear modulus [i: 



^bT , 2 



V B 



-{e 



1 



201 
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This relation depends also on the ratio Vb/V, as the 
energy required for the rotation of a disk, which is not 
embedded in a surrounding medium, tends to zero. The 
analysis of the strain variable 20 offers thus an indepen- 
dent, direct route to the determination of the shear mod- 
ulus in an embedded system. 

These considerations show that in order to obtain ac- 
curate elastic moduli from the analysis of the strain fluc- 
tuations the relative size of the analyzed system to the 
complete, finite system should to be known. Neverthe- 
less for the case of the colloidal crystal [ll[ the situation 
is close to the limiting case of ^ — » 0. Here the influ- 
ence of the surrounding medium on the analyzed system 
is dominant. The strain variables e\ and e2e can be used 
to extract the elastic moduli. Thus the two-dimensional 
colloidal crystal, as discussed in detail in [lj, [l8j is an 
example for a completely embedded system. In contrast 
to this in simulations the complete system can be ana- 
lyzed, which corresponds to the limiting case of ^y- — > 1. 
For this case each of the strain-strain correlation func- 
tions of the strain variables e±, e2 and e% give directly 
access to the corresponding elastic moduli dj. The effect 
of the embedding medium can be visualized by looking 
at a statistical sum rule, as will be discussed in the next 
paragraph. 



3. Analysis of a statistical sum rule 

The sum rule for the generalized susceptibility provides 
another way of extracting the elastic moduli from the 
strain-strain correlation functions. The coarse-grained 
system represents a homogeneous continuum and is thus 
translationally invariant. For such systems the suscepti- 
bilities are directly related to the correlation functions. 
In Fourier space this reads xr(fc) = /3G(fc). From this 
the static susceptibility sum rule follows [17| : 



XT =lim XT(k)=pG(k)\^ g = p 

k^O 



dr G(r) 



Thus an integration of the correlation functions in real 
space yields directly the elements of the compliance ten- 
sor Sijki , which correspond to the generalized susceptibil- 
ities. The compliance tensor is the inverse of the stiffness 
tensor Bijki. In the case that no external stresses act on 
the system this is equivalent to the tensor of the elastic 
constants Cijki [12] ■ The S^ki obtained from such an 
analysis depend on the integration volume Vb- Thus in 
order to obtain systems-size independent values an ad- 
ditional finite-size scaling analysis should be employed. 
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FIG. 9: (Color online) Compliances Su as they are obtained 
from the sum rule as a function of the ratio of analyzed volume 
Vb = L% to the complete simulation volume V — L 2 . Shown 
are the results of Monte Carlo simulations of a harmonic tri- 
angular lattice with N = 5822 particles and spring constant 
/3a 2 / = 200/^3. Lines are a fit to the data with the formulae 
given in the text, while the dotted lines are a guide for the 
eyes, a) NpT ensemble with periodic boundary conditions, b) 
NVT ensemble with periodic boundary conditions, c) NVT 
ensemble with open boundary conditions. Here dashed lines 
show the value for Lb/L = 1 as a comparison. 



Figure [5] shows the compliances Su (i = 1, 2, 3 and 20) 
as a function of the ratio of the integration volume Vb to 
the complete simulation volume V, i.e. Vb/V = Lb/L, 
as they are obtained from the simulation data of the har- 
monic triangular crystal at zero external pressure in the 
NpT and NVT ensemble with different boundary con- 
ditions. A comparison shows directly how the choice 



of ensemble and the choice of boundary conditions in- 
fluences the results. These are so called explicit and 
implicit finite-size effects (23|. In addition this analy- 
sis visualizes the effects of the embedding medium on 
the compliances Su . Figure [5] a) shows the results from 
simulations in the NpT ensemble with periodic bound- 
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ary conditions. The complete system contains N = 5822 
particles, that are connected via springs of spring con- 
stant (3a 2 f — 200/V3- The strain-strain correlation func- 
tions are directly related to the elastic moduli in this 
ensemble, due to the fact that the volume itself fluctu- 
ates. Thus one can obtain the elastic moduli directly 
from the Sn at Lb/L = 1. For the considered system 
we find: s[[ ] = 0.01025, S$ = 0.02050, S$ = 0.00512 

and S^glg —> 0, resulting in a\ — 97.5, a2 = 48.7 and 
as = 195.2. These values lie within 2.5% of the expected 
values. The accuracy of this approach compares to that 
of the methods for the calculation of the elastic moduli 
discussed before. From the considerations in IIII D II one 
expects the following functional dependence of Sn on 
L B /L: 

s[i b) = V B {e\)/{k B T)= [^ + M (1-(L|/L 2 ))]- 1 

The black solid line in figure [pj a) is a fit with this equa- 
tion to the data (crosses). From the fit parameters the 
following elastic moduli are extracted: a\ = K = 93.4 
and n = da = 45.8. Figure [5J a) shows clearly the in- 
creasing impact the surrounding medium has on Sn as 
Lb/L diminishes. In the limit Vb/V — > it yields the 
sum of the elastic moduli K and fi. It is apparent from 
figure 03 a), that as soon as Vb/V < 1, the compliances 
S22 and S33 cannot be directly related to the shear mod- 
ulus any more. 

The considerations in IIII D 2l suggest, that the compli- 
ance S2020 should diverge as Vb/V — ► 1. This relates to 
the fact, that the energy needed for rotating an embedded 
disk goes to zero as the embedding material is removed. 
This divergence cannot be seen in the simulation data, 
as in simulations with periodic boundary condition the 
system as a whole cannot rotate. The fact, that there is 
no divergence of S2020 for Vb/V — > 1, is therefore an im- 
plicit finite-size effect. In order to extract the shear mod- 
ulus from the compliance 82026 a polynomial in (Lb/L) 
was fitted to the data (open triangles). From the limit 
Vb /V — ► the shear modulus is extracted: fj, = 48.9. 

In the NVT ensemble with periodic boundary condi- 
tions the compliances as a function of Lb/L exhibit a 
different dependence on Lb/L, as figure 0b) shows. An 
unstrained state of the triangular lattice is analyzed in 
these simulations, therefore the integral of the correla- 
tion functions over the complete system goes to zero and 
gives no access to the elastic moduli. This is an explicite 
finite-size effect. Nevertheless from the limit V B /V — ► 
one can extract the elastic moduli as in the case of the 
simulations in the NpT ensemble with periodic boundary 
conditions from the compliances Sn (crosses) and S2926 
(open triangles). Fits with a polynomial in (Lb/L) are 
plotted as solid lines in figure[5]b). From the case of max- 
imum embedding one extracts a\ + 02 — K + fi = 148.3 
from Sn and ci2 = 51.5 from 82020 in figure[H]b). These 
values compare to the values obtained by different meth- 
ods as they are given in table [ITJ 

The compliances S^ B shown in figure [5] c) are ob- 
tained from data of simulations in the NVT ensemble 



with open boundary conditions as they were presented 
in [TH |. These show in contrast no systematic depen- 
dence on Vb /V. The maximum analyzed volume, which 
will for this case be denoted by F = L 2 , is approxi- 
mately one fourth of the complete system volume. Av- 
eraging over the positions of origin, as it is done in the 
calculation of the strain-strain correlation functions in 
the system with open boundaries, results in an averag- 
ing over sub-systems with partial to complete embed- 
ding. For this type of averaging the k = values of all 
considered strain-strain correlation functions give access 
to the elastic moduli of the system [ll[. The effect of 
this type of averaging shows up most prominently in the 
fact that S2028 does not tend to zero for Lb/L — > 1, 
but approaches the value of 522- Extracting the elas- 
tic moduli from the compliances for Lb/L = 1 in fig- 
ure M c) yields s{^ = 0.00940 ->■ ax = 106.9, S% } = 
0.02100 -> a 2 = 47.6, S$ = 0.00525 -> a 3 = 190.4 and 
S 2ele = 0.02160 -> a 2 = 46.4. The accuracy of these 
values is the same as in [ll[ . The deviation from the the- 
oretical values is larger in this the finite system 
with open boundary conditions is influenced in its elastic 
properties by the missing, stabilizing bonds for particles 
at the surfaces. 



IV. THE STATISTICS OF NON-AFFINE 
FLUCTUATIONS 

The coarse-graining process described in section IIIII 
projects the configurations generated by our microscopic 
Hamiltonian onto strain fields which are smooth over dis- 
tances larger than the coarse-graining length A. It also 
generates a conjugate noise [15| which represents those 
fluctuations which cannot be captured during coarse- 
graining. This is easily understood once it is realized 
that, coarse-graining retains only that part of the par- 
ticle displacements u = f — R in a configuration which 
can be obtained from the reference lattice R by an affine 
transformation: r = (1 + e)R. An affine transformation 
constrains all parallel lines in the reference lattice to re- 
main parallel, which is clearly impossible to satisfy for 
an arbitrary configuration coarse-grained over volumes 
larger than an unit cell. Indeed, the quantity x as de- 
fined in Eq. ([7]) has the dimension of Length 2 and scales 
as A 2 . This may be seen by comparing figure [10] (a) 
and (b) . Figure [TU] shows the probability distribution 
P(x) and its scaling behavior for various choices of pa- 
rameters. While [TU] (a) shows a clear dependence of the 
amount of non-affinity on the coarse-graining length A, 
figure [TD1 (b) shows a collapse of the distributions for the 
scaled quantity x/A 2 for A > 2.2. This corresponds to 
a minimum of 18 neighbors to the central particle, that 
are taken into account in the calculation of the strain 
field via the minimization of \- These distributions show 
a constant offset from x/A 2 = 0. By contrast A = 1.3 
shows no such offset, meaning that for the calculation of 



14 



the affine strain field only the minimal neighborhood, i.e. 
the 6 nearest neighbors, allows for a global minimization 
of x- In addition the probability distributions of the non- 
affine parameter x a l s ° scale with with the spring con- 
stant / (figure [10] (c)) and in simulations run at various 
hydrostatic, external pressure p with the resulting aver- 
age density (g) (figure [TDl (d)). One can therefore obtain 
the probability distribution for x for anv inverse temper- 
ature /3, spring constant /, density p and coarse-graining 
length A from a generalized extreme value probability 
distribution function. This master curve is P{X) — 



with 



(1 + 0.27 • z y-^-«>-L e -(i+o.27.*) j /2 012 

z = (X — 3.127) /2.012, which we obtain from a fit to the 
simulation data, with the scaling variable X = x/3///?A 2 , 
independent of system size N and the choice of ensemble. 

We show below that features of -P(x), like the de- 
pendence on the spring constant /, may be rationalized 
within a simple "cell model" calculation. In this model 
each particle is assumed to fluctuate within the cage of 
its 6 nearest neighbors which suffers, at most, an affine 
distortion (see figure [TT|). The only source of non- affinity 
comes from the displacement of the central particle from 
its equilibrium position. For such a subset of configu- 
rations, one may simply decompose each configuration 
{r} as that obtained by an affine transformation plus a 
non-affine displacement s of the central particle within an 
undistorted hexagonal cell. The non-affinity parameter 
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FIG. 10: (Color online) Plots of the probability distribution 
of the non-affinity parameter x from simulations in the NpT 
ensemble with periodic boundary conditions and N = 3120 
particles, (a) Plot of the probability distribution of P(x) vs x 
for various coarse- graining length A. (b) The probability dis- 
tribution of x/A shows a data collapse for A > 2.2. (c) Data 
for various spring constants / collapse onto each other. The 
prediction from a simple cell model (black line) is shown for 
comparison, (d) Data from simulations at various hydrostatic 
pressure p scale with the resulting average density (g). 
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FIG. 11: (Color online) Example configuration (a) consisting 
of a central particle and neighboring particles i — 1,6 which 
can be decomposed into a purely affine deviatoric distortion 
(b) together with a non-affine displacement s of the central 
particle (c). 



X may then be calculated to be, 

6 2 2 



m— 1 i— 1 
6 2 



= EE( r ™- r S"^-^)) = 6s2 

Within this approximation the ensemble average con- 
tributes to the Lindcmann parameter I, as (|it| 2 ) = 
<(C« nf + s) 2 > = (ul ffine ) + { X )/G + 2{K mns \*\s\) = IV. 
The Lindemann ratio depends on the stiffness of the solid 
and grows as the melting point is approached. Within 
this model the energy of these configurations can be cal- 
culated to be, 



E 



f 



ARi 



f 



-2s • ARi - 2Ai? 2 y / l + s 2 /Ai? 2 - 2 ARi ■ s/Ai? 2 j 
« 3/s 2 

Here we used an approximation of the square root up to 
0(4) and the abbreviation ARi = R4 — Rq. Thus the en- 
ergy related to the non-affinity x of the central particle is 
Eceii = E/3 = /x/6. With this energy contribution it is 
straight forward to calculate the probability distribution 
of X, 

P[X) = C J dse-^ s2 S(x'(s)~ X ) 
2 (Pf 



(7) 



where <f> = @<fx an d C, the normalization constant. In 
figure [10] (b) we have plotted P(x)/(/3//6) from this cell 
model together the scaled distributions obtained from 
our simulations. The data collapse of the distributions 
from simulations with various spring constants is in ac- 
cord with the scaling in /3//6 as expected from the simple 
cell model. The details of the shape of the distribution 
function can not be captured completely. As is to be ex- 
pected in this simple model, the contributions of large x 
are slightly overestimated. 

How docs the presence of x influence strain correla- 
tions? To see this we assume that, at least for small x 
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FIG. 12: (Color online) a) The autocorrelation function of the 
non-affinity parameter \ f° r the harmonic triangular lattice 
with spring constant /3a 2 / = 200/\/3 in the NpT ensemble 
with periodic boundary conditions along the x and the y axes. 
The red line is a fit to an exponential form, b) A surface plot 
of its Fourier transform G xx (k). 



the total strain obtained by fitting an arbitrary configu- 
ration to an affine transformation contains an affine part 
which would have been the only result if \ were zero, and 
a x dependent non-affine part which may be expanded 
as a series in powers of x, namely, 

<%•(%) = £ij+^2t P X p 
v 

This decomposition is more general than what is sug- 
gested above, and it is customary, in theories of solid plas- 
ticity to decompose the total strain into elastic (affine) 
and plastic (non-affine) parts efj = e°» + [24| . To low- 
est order in \ therefore, 

(e(x;0)ye( X ;r) w ) = (e°(0)ye°(r)«) + f?<*(0)x(f)) 

where we have used the fact that the coarse-graining pro- 
cess projects the displacements into mutually orthogonal 
subsets [HI, [26[ so that one can ignore all correlations 
between A and x- I n figure fT2l a) we have plotted cuts 
showing the decay of G xx (r) = (x(O)x(r')) along the x- 
and y- axis. The function G xx is isotropic and decays 
rapidly to zero over a length scale comparable to A. This 
suggests that x behaves as a "delta" correlated white 
noise with a probability distribution given by Eq. [JJ 
Again, this is consistent with our identification of x with 
the Lindemann ratio, the microscopic, random, thermal 
fluctuations of individual particles are, indeed, expected 
to be uncorrelated with each other. Given the form of 
G xx one expects such fluctuations to contribute only a 
background term (compare figure [12] b)) to the strain 
correlations in Fourier space. 

We have shown in this section that the coarse-graining 
process by which affine strains may be extracted from 
microscopic particle configurations also generates a ran- 
dom white noise consisting of non-affine particle displace- 
ments. For a harmonic solid this is simply related to 



the Lindemann parameter, which, in turn, depends ulti- 
mately on the strength of the interactions. 

What is the general implication of this to the study of 
elasticity and rheology of complex solids? The current 
picture of the mechanism of relaxation in amorphous 
materials indicates that there are two main competing 
processes involved [25j. Over small time scales the 
system fluctuates within local minima in the free energy 
landscape making transition between such basins of 
attraction over longer time scales. We have shown here 
that harmonic fluctuations within local minima gener- 
ates a well characterized contribution to the non-affine 
displacement Xj therefore any "extra" contribution to 
X arises exclusively from these inter-basin transitions. 
Thus our analysis may be used as a tool to distinguish 
between these two kinds of relaxations in complex solids. 



V. CONCLUSIONS 

We have shown in this paper how the analysis of parti- 
cle configurations of two-dimensional soft solids gives ac- 
cess to a wealth of information on the local and non-local 
elastic properties. Since the harmonic solid analyzed here 
is the most generic conceivable our work has the poten- 
tial to serve as a template for further research in this 
direction. The properties of the strain-strain correlation 
functions have been discussed in great detail and various 
methods of how to extract the elastic moduli from their 
analysis were presented. Furthermore we determined and 
discussed the effects of external pressure and an embed- 
ding medium, the proper treatment of which is essential 
for experimentalists seeking to use our methods for an- 
alyzing mechanical behavior of soft matter. The impli- 
cations of our work particularly for the understanding of 
non-affineness in solids is significant, because our study 
allows one to classify non-affine fluctuations in any sys- 
tem into "trivial" (in the sense of being present even in an 
ideal harmonic solid) and non-trivial components. In the 
future, we shall use these procedures to study metasta- 
bility in solids undergoing phase transitions and plastic 
behavior of solids under large external stresses. 



Acknowledgments 

We acknowledge useful discussions with K. Binder, R. 
Messina and M. Rao. This work was funded by the 
Deutsche Forschungsgesellschaft (SFB TR6/C4). Grant- 
ing of computer time from HLRS, NIC and SSP is grate- 
fully acknowledged. One of us (SS) thanks the DST, 
Govt, of India for support. 



16 



[1] A. Ycthiraj, A. van Blaaderen, Nature, 421,513 (2003). 
[2] P. Habdas, E. R. Weeks, Curr. Opin. Colloid Interface 

Sci., 7, 196 (2002). 
[3] K.J. Strandburg, Rev. Mod. Phys., 60, 161 (1988); K.J 

Strandburg, ibid., 61, 747 (1989). 
[4] S. Sengupta, P. Nielaba and K. Binder, Phys. Rev. E, 

61, 6294 (2000). 
[5] K. Binder, S. Sengupta and P. Nielaba, J. Phys.: Con- 
dens. Matter, 14, 2323 (2002). 
[6] H.H von Griinberg, P. Keim, K. Zahn and G. Maret, 

Phys. Rev. Lett., 93, 255703 (2004). 
[7] A. Wille, F. Valmont, K. Zahn and G. Maret, Europhys. 

Lett., 57, 219 (2002). 
[8] S. Sengupta, P. Nielaba, M. Rao and K. Binder, Phys. 

Rev. E, 61, 1072 (2000). 
[9] K. Zahn, A. Wille, G. Maret, S. Sengupta and P. Nielaba, 
Phys. Rev. Lett., 90, 155506 (2003). 
[10] R. Maranganti, P. Sharma, Phys. Rev. Lett., 98, 195504 
(2007); R. Maranganti, P. Sharma, Journal of the Me- 
chanics and Physics of Solids, 55, 1823 (2007). 
[11] K. Franzrahe, P. Keim, G. Maret, P. Nielaba and S. Sen- 
gupta, Phys. Rev. E, 78, 026106 (2008). 
[12] M.L. Falk, J.S. Langer, Phys. Rev. E, 57, 7192 (1998). 
[13] A. Lemaitre, Phys. Rev. Lett., 89, 195503 (2002). 
[14] C.E. Maloney, M.O. Robbins, J. Phys.: Condens. Matter, 

20, 244128 (2008). 
[15] P.M. Chaikin, T.C. Lubensky, Principles of condensed 



matter physics (Cambridge University Press, Cambridge, 
UK, 1995). 

[16] D.G.B. Edelen in Continuum Physics 4, A.C. Eringen, 
eds, (Academic Press, New York, 1976). 

[17] N. Goldenfeld, Lectures on Phase Transitions and the 
Renormalization Group (Westview Press, Boulder, Col- 
orado, USA, 1992). 

[18] K. Franzrahe, P. Nielaba, A. Ricci, K. Binder, S. Sen- 
gupta, P. Keim and G. Maret, J. Phys.: Condens. Mat- 
ter, 20, 404218 (2008). 

[19] P. Keim, G. Maret, U. Herz and H.H. von Griinberg, 
Phys. Rev. Lett., 92, 215504 (2004). 

[20] M. Parrinello, A. Rahman, J. Chem. Phys., 76, 2662 
(1982). 

[21] D.R. Squire, A.C. Holt and W.G. Hoover, Physica, 42, 
388 (1969). 

[22] D.C. Wallace, Thermodynamics of Crystals (Dover Pub- 
lications Inc., Mineola, NY, 1998). 

[23] F. L. Roman, J. A. White and S. Velasco, J. Chem. 
Phys., 107, 4635 (1997); F. L. Roman, J. A. White, A. 
Gonzalez, S. Velasco J. Chem. Phys., 110, 9821 (1999). 

[24] J. Lubliner, Plasticity Theory (Dover Publications Inc., 
Mineola, NY, 2008). 

[25] K. L. Ngai, G. B. Wright, eds, Relaxations in Complex 
Systems (NRL, Washington, DC, 1985). 

[26] H. Mori, Progress of Theoretical Physics, 33, 423 (1965). 



